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Abstract. We have developed a fast numerical 2-D model of galaxy disk 
evolution (resolved along the galaxy radius and azimuth) by adopting a scheme 
of parameterized stochastic self-propagating star formation. We explore the 
parameter space of the model and demonstrate its capability to reproduce 1-D 
radial prohles of the galaxy M33: gas surface density, surface brightness in the 
i and GALEX FUV passbands, and metallicity. 
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1. INTRODUCTION 

Numerous recent resolved photometric and spectroscopic surveys of galaxies 
allow us to explore a two-dimensional (2-D) structure of galaxy disks (resolved 
along the galaxy radius and azimuth) in detail. Studies of disk structures have 
proved them to be feature-rich (e.g., Bastian et al. 2007; Gieles et al. 2008; 
Bastian et al. 2009; Sanchez et al. 2010). On the other hand, simulations of galaxy 
disks with hydrodynamical models are computationally costly and dependent on 
the parametrization of sub-grid physics and even on the implementation methods 
(Scannapieco et al. 2012). Additionally, in the cases of late-type spirals hne tuning 
of the models is needed (Roskar et al. 2014, and references therein). An increased 
spatial resolution of the simulations could be the key to the problem, however, 
this would make galaxy models computationally even more daunting (Guedes et 
al. 2011). 

In the past, there was an attempt to create fast semi-analytic 2-D models 
of galactic disks, based on a stochastic self-propagating star formation (SPSF) 
scenario (Seiden & Gerola 1982, and references therein). The SPSF approach was 
shown to be suitable to form flocculent spiral arms reminiscent of late-type spirals 
(Gerola & Seiden 1978). The properties of dwarf galaxies (Gerola et al. 1980), 
as well as larger spirals (Feitzinger et al. 1981), were explained by employing 
the SPSF scenario. Despite further developments of the grid-free SPSF variants 
(Jungiewert & Palous 1994; Sleath & Alexander 1995), galaxy modeling shifted 
towards N-body and hydrodynamical approaches. 

In this study we present an updated version of the 2-D stochastic galaxy disk 
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model by Mineikis & Vansevicius (2010). The most important new features of the 
updated version are the implementation of the radial profile of gas accretion onto 
the disk, based on the present-day mass distribution in the galaxy M33, and the 
stellar “dispersion” across the disk. These features allow us to produce realistic 
radial profiles of the surface brightness. We also improved the model by minimizing 
a number of free parameters and designing it for fast generation of the results of 
resolved stellar photometry. 

In order to verify the model, we explored an extensive grid of parameters by 
comparing the model generated radial profiles of one-dimensional (1-D) galaxy 
disk with the observed ones. For the study we used the well explored Local Group 
galaxy M33. The relatively unperturbed stellar disk (Ferguson et al. 2007) and its 
favorable inclination together with the small distance - 840 kpc (Freedman et al. 
2001), as well as the low Galactic extinction in its direction, makes this galaxy one 
of the best laboratories for studies of the SPSF scenario. M33 has the characteristic 
flocculent spiral arms inherent to the late-type disk galaxies. Although two spiral 
arms are visible in the near infrared (Regan & Vogel 1994), they do not dominate 
the present-day star formation pattern, as can be seen on the far ultraviolet and 
Ha images (Thilker et al. 2005). The abundance of observational data makes the 
M33 galaxy the best target to test our 2-D model. 

The article is organized as follows. In Section 2 we describe the model, in 
Section 3 the model parameters are calibrated for the galaxy M33, a discussion of 
the simulation results is given in Section 4 and, finally, in Section 5 the conclusions 
are presented. 

2. THE MODEL 

2.1. Disk geometry and time step 

The disk model is divided into Vr rings of equal width (Mineikis & Vansevicius 
2010). Each ring is subdivided into cells, which are the smallest elements of the 
model. This subdivision follows the rule that the ring with the running number 
i has 6 • i cells, producing a total of — 1) -I- 1 cells in the disk. This 

subdivision results in equal area cells, except for the central one, which is smaller 
by a factor of 3/4. The physical size of the disk model is defined by the physical 
size of cells and the number of rings. 

We assume the size of cells to be comparable to the typical size of OB associ¬ 
ations. According to the recent findings by Bastian et al. (2007), star formation 
(SF) is a hierarchical scale-free process, highly dependent on the definitions. Nev¬ 
ertheless, our model is not very sensitive to the adopted cell size, therefore, in this 
study, a cell size is set to dc = 100 pc. 

The model integration time step (Afj) corresponds to the SF propagation time 
across the cell. We assume that the SF propagation velocity is usf = 10 km/s 
(Feitzinger et al. 1981), i.e., it corresponds to the typical speed of sound in the 
interstellar matter within a cell, therefore: 


dc ^ / VSF \ 
100 pc \10 km/s / 


Ati = 10 Myr • 


-1 


( 1 ) 
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2.2. Gas accretion 

Although mass accretion histories of dark matter (DM) halos are highly vari¬ 
able (McBride et al. 2009), they tend, on average, to obey simple relations. 
McBride et al. (2009) studied the Millennium simulation data and found a simple 
empirical fit to the average DM accretion rate. The fit relates the mean accretion 
rate of DM to a given halo mass at a given redshift. Using the higher resolution 
Millennium II simulation, Fakhouri et al. (2010) confirmed the validity of the fit 
and extended it to smaller halo masses. We use the results by Fakhouri et al. 
(2010) to prescribe the DM and, correspondingly, baryonic mass (BM) build-up of 
the model. 

The gas falling into the DM halo is distributed in the thin disk. The radial 
profile of the accretion is assumed to be a scaled version of the total BM radial 
profile of the present-day galaxy disk: 


Ag(d t ) 


B{r) 

2tt B{r)rdr 


• Adm)^) • P, 


( 2 ) 


where AG{r,t) is the radial profile of gas accretion rate, B{r) is the radial profile 
of BM density in the galaxy disk, ADM(i) is the rate of DM accretion on the 
galaxy halo, and /3 = Dbm/Udm is the primordial ratio of BM to DM. Such a 
definition of accretion guarantees the build-up of the present-day disk over the 
galaxy’s life-time. The metallicity of the accreted gas is set to Za = 0.0001. 


2.3. Star formation 

The star formation process in the disk is modeled by discrete SF events occur¬ 
ring stochastically in the cells. The cell can experience a SF event spontaneously 
(probability Ps) or be triggered (probability Pt)- 

• Spontaneous SF probability Pg: this parameter defines SF events occurring 
stochastically in the cells without any external influence. The spontaneous SF 
sustains SF activity in the disk. In this study we set Ps to be very small, i.e., ~I 
SF event per model integration time-step Ati over the whole disk. Additionally, 
we assume that spontaneous SF events are more probable in the cells possessing 
higher gas density: Ps c< Eq, where Eg is a surface gas density within a cell. 

• Triggered SF probability Pp: this parameter represents a chain of complex 
processes in the cell’s interstellar medium (molecular cloud) after the SF event 
has occurred. A molecular cloud will undergo imminent disruption by an energetic 
feedback of stellar winds, expanding H II zones and supernovae explosions. Despite 
a negative SF feedback locally (on a scale of cell), SF could be triggered on a larger 
scale (i.e., in neighboring cells). For the cell i, which experiences a SF event in 
a time step tk, the neighboring cell j, being in direct contact with the cell i, can 
experience a SF event during the next time step tk+i with a probability of Ptj- 

At the start of galaxy simulation the SF in the disk is inhibited by the critical 
gas surface density Ec. If the gas surface density in a cell is below critical, the 
probability of SF event is reduced by a factor of Eg/Ec. As the simulation evolves, 
the gas density in the disk increases and SF becomes more probable starting from 
the inner disk parts. In the outer disk. Eg remains small, compared to Eg, for a 
longer time and defines the extent of the SF disk. We note that Eg defines the 
local density, i.e., not averaged azimuthally, as this parameter is usually derived 
from observations. 
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Table 1. PEGASE-HR parameters used to generate SSPs. 


Parameter 

Value 

Reference 

Stellar library 

low-resolution 

Le Borgne et al. (2004) 

Initial mass function 

corrected for binnaries 

Kroupa (2002) 

Fraction of close binaries 

0.05 

default PEGASE-HR value 

Ejecta of massive stars 

type B 

Woosley & Weaver (1995) 

Nebular emission 

true 

PEGASE-HR value 


During a SF event, a fraction of gas available in a cell is converted to stars. 
This fraction is called the SF efficiency, SEE: 


"'(lOMo/pcO ’ 

where e and a are free parameters. In order to avoid generating unrealistic stellar 
populations of extremely low mass or SEE exceeding 100%, we set the lower and 
upper SEE cut-offs to 0.05% and 50%, respectively. Additionally, if the stel¬ 
lar population, born in a cell, is less massive than 100M©, we assume that the 
starburst does not trigger neighboring cells because the population is, on average, 
lacking supernovae. 

2.4- The evolution of cells 

Each stellar population, formed in a cell during a SF event, could be well 
represented by a simple stellar population (SSP). To track the entire SF history of 
a cell, the value of mass of each particular stellar population is stored in a 2-D (age 
vs. metallicity) array, S. The step in age, tg, can be set in accordance to the size 
of a cell (10 Myr throughout the paper). The step in metallicity (0.1-0.3 dex) is 
predefined by the used dataset of interpolated isochrones. To track the evolution 
of SSPs self-consistently, we employ the package PEGASE-HR (Le Borgne et al. 
2004; see Table 1 for the parameters used). 

The stars from any SSP, formed in a cell, can move to the neighboring cell, 
therefore, arrays Sj of the neighboring cells have to be modified accordingly. Stellar 
spread to neighboring cells is implemented as a “dispersion” process. The evolution 
of the stellar content in the cell Si is given by 


3 

where dti is the mass of a newly formed SSP, Ds is the “dispersion” constant, 
is the length of borderline between cell i and its neighbors j. To prevent star flow 
anisotropies near the center of galaxy, which occur due to change in the length 
ratios of the cell sides, we correct corresponding A to keep star flows through each 
side of the cell equal. 

The evolution of gas content in the cell Gi is given taking into account SSP 
evolution: 

AG- 

= + E ° (5) 

j k I 


where Ai is the accreted gas mass, Fij are the gas mass flows (see Section 2.5) 
between cell i and neighboring cells j, and R is the array representing the mass 
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fraction of stellar populations returned to gaseous content, i.e., the last terirQ 
represents the total gas mass returned to gas pool of the cell i by stellar populations 
evolving within this cell (the indices k and I denote the age and metallicity). 

The metallicity (Zi) evolution of the gas content in the *-th cell is given by 


=-^^-z,+a-za-J 2 ■ Zj +° ( 6 ) 

i fc I 

where Za represents the metallicity of the accreting gas, j denotes the neighboring 
cells, and Z is the array representing yields of metals returned to the gas pool of 
the i-th cell by stellar populations evolving within this cell. 

2.5. Gas flows 

Due to the small size (~ 100 x 100 pc) of cells covering the galaxy disk, gas 
masses are expected to flow beyond the cells in one time step. The main structures 
causing gas movement are superbubbles inflated by SF feedback. Following the 
formulation by Castor et al. (1975), on time-scale of 10 Myr superbubbles reach 
~ 100 pc in size for a typical gas density of 20 Mq/pc^ and SFE values of a few 
percent. Therefore, a cell experiencing a SF event on a time-scale of 10 Myr will 
be filled with hot tenuous gas, causing a large part of the gas to move out of the 
cell. 

After 40 Myr, when supernovae vanish and the superbubble inflated cavity 
cools down, the gas returns to the cell. The force driving the gas to refill the cell 
is the random motion of HI gas clouds having typical velocities of ~ 10 km/s (see, 
e.g., Corbelli & Schneider (1997) for the galaxy M33). Based on the arguments by 
Hunter & Gallagher (1990), we estimate the refill time-scale t = dc/{‘2-10 km/s) ~ 
50 Myr. 

We implemented both types of gas flows in the galaxy disk model: expulsion 
by a SF event and refilling. 

• Gas expulsion: this is implemented simply by moving gas from the cell i, 
experiencing a SF event, to the neighboring cells j. The gas flows only to those 
neighbors which have had no a SF event for the past 40 Myr: 


F- — \ 

-^ 2,7 — 


0 ^SF J < 40 Myr 
Gi Asf J > 40 Myr, 


( 7 ) 


where gas flows between cells depend on the length, of the borderline, 
and Gi is the gas mass remaining in the cell i after the SF event. The flow occurs 
within one time step after the SF event. 

• Gas refilling: this is implemented under the assumption that equilibrium 
gas distribution in the galaxy follows the present-day baryonic matter distribution 
B{r), i.e., any deviation from it, {G{ri)/G[rj) B{ri)/B{rj)), generates gas flows 
between cells on the time-scale r: 


F- - — \ - 


G{r,) ■ Bin) - G(r,) ■ B{r,) 
B{r,) + B{rj) 


( 8 ) 


^Also known as Hadamard (or Schur) product, for two matrices of the same dimension, 
{B o C)k,i = • Ck,i 
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Table 2. The parameters of the 2-D galaxy model. 


Parameter 

Notation 

Value 

Fixed parameters 

Age 

- 

13 Gyr 

Baryonic mass 

- 

1.1 • 10^° Mq 

Disk radius 

- 

12 kpc 

Disk rotation 

- 

Corbelli & Salucci (2007) 

Cell size 

dc 

100 pc 

SF propagation speed 

VSF 

10 km/s 

Accretion rate 

ADM(i) 

Fakhouri et al. (2010) 

Critical gas surface density 

Sc 

8 Mq/ pc^ 

Gas refilling time-scale 

T 

50 Myr 

Star “dispersion” constant 

Ds 

200 pc^/Myr 

Spontaneous SF probability 

Ps 

1 SF region per Atj 

Varied parameters 

Triggered SF probability 

SF efficiency 

SF power index 

Pt 

e 

a 

0.28-0.44 
(0.018-5.7)-10-2 

1.5-2.5 


3. THE MODEL SETUP EOR M33 

The model parameters used for the M33 galaxy simulation are given in Table 2. 
The model disk (radius of 12 kpc) comfortably engulfs the star forming disk of M33. 
The presence of the RR Lyr variables in M33 (Pritzl et al. 2011 and references 
therein) implies its old age, therefore, the galaxy formation was assumed to start 
13 Gyr ago. The critical gas density for SF (Sc) was suggested to be in the range 
3-10 Mc/pc^ (Schaye 2004). We have found that lower end values do not change 
the simulation results significantly and adopted Sc = SMc/pc^. The 1-D radial 
profiles are not very sensitive to the stellar “dispersion” parameter, therefore, it 
was set to 200 pc^/Myr. The effects of this “dispersion” parameter will be analyzed 
in the subsequent paper devoted to the study of properties of 2-D galaxy models. 

An assumption was made that the radial profile of gas accretion is a scaled 
version of the radial profile of total baryonic mass (gaseous and stellar) in the 
galaxy. To avoid computationally expensive iterative fitting of baryonic mass 
for each model, we derived the radial profile by performing the decomposition of 
the galaxy’s rotation curve. The observational data of the rotation curve were 
taken from Corbelli & Salucci (2007). The rotation curve was decomposed into 
four mass components: gaseous disk, stellar nucleus, exponential stellar disk and 
pseudo-isothermal DM halo. 

The contribution of the gas disk to the rotation curve was calculated by co¬ 
adding HI and H 2 data from Corbelli & Salucci (2000) and Heyer et al. (2004), 
respectively. Helium mass fraction was taken into account by multiplying hydrogen 
gas mass by a factor of 1.33. Corbelli (2003) has established that the mass of the 
M33 stellar nucleus is in the interval (0.3 — 8) • 10® Mq, therefore, in order to narrow 
a free-parameter space, we assume that the stellar nucleus is a point source with 
a mass of 10® Mq. To avoid the effects of the extended nucleus, we did not fit the 
innermost 0.5 kpc part of the rotation curve. 
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R, kpc 

Fig. 1. Circular speed versus radius for the galaxy M33. The black dots with error 
bars indicate the observed rotation curve (Corbelli & Salucci 2007). The red line is the 
best-htting model rotation curve (see the text for the model description). Individual 
contributions from the components are also shown: pseudo-isothermal DM halo (cyan 
squares), exponential stellar disk (blue stars), gaseous disk (green diamonds) and stellar 
nucleus (magenta triangles). 

We fitted the following two components: 

• a pseudo-isothermal DM halo with the core radius Rc and the asymptotic 
rotation velocity Vii{oo), e.g., de Blok et al. (2008); 


• an exponential stellar disk with the central surface density Eg and the scale 
length Re. The stellar disk was assumed to have a scale height of 100 pc. 

Dynamical contributions of gaseous and stellar disk components were computed 
using the task rotmod within the software package GIPSY (van der Hulst et al. 
1992) which implements the method of Cesartano (1983). For each pair of the 
parameters describing the stellar disk, Eg & ’^6 minimized by fitting the 

DM radial profile. 

The best fitted rotation curve with Re = 1.57kpc, Eg = 470 Mq/pc^, Rc = 
10.2kpc, and ^^(oo) = 184km/s is shown in Fig. 1. The derived scale length of 
stellar disk is in agreement with the values derived photometrically, 1.4 kpc in the 
7f-passband for the inner 4 kpc of M33 (Regan & Vogel 1994) and 1.5-1.6 kpc in 
the Spitzer passbands at 3.4 and 4.5 (Verley et al. 2009). The best-fitted 
stellar and gaseous disks were used to define the radial gas accretion profile of 
M33. 

The grid of models was computed by varying the main SF parameters: e, a 
and Ft- The grid is composed of 17 linear steps in Px, 5 linear steps in a and 
11 logarithmic steps in e. In order to reduce stochastic effects, every model is 
presented as an average of six snapshots spaced by 200 Myr during the last I Gyr. 
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R. kpc R. kpc 

Fig. 2. Comparison of the observational data on the galaxy M33 with the models. 
The gas surface density (panel a) is computed by co-adding the radial profiles of HI 
(Corbelli & Salucci 2000) and H 2 (Heyer et al. 2004) gas surface density. The surface 
brightness radial profile in the i passband (panel b) (Ferguson et. al. 2007) is corrected 
for internal extinction using the radial extinction profile by Munoz-Mateos et al. (2007) 
and assuming the LMC-like extinction law (Gordon et al. 2003). The surface brightness 
radial profile in GALEX FUV passband (panel c) (Munoz-Mateos et al. 2007) is also 
corrected for internal extinction. The metallicity measurements of blue supergiant stars 
(panel d) are taken from Urbaneja et al. (2005) and U et al. (2009). All radial profiles of 
observed surface brightness are de-projected adopting 54° for the galaxy disk inclination. 


4. RESULTS 

Comparisons of the models with observational data for the galaxy M33 are 
presented in Fig. 2. The models are compared with 1-D radial profiles of: gas 
surface density (HI from Corbelli & Salucci 2000, H 2 from Heyer et al. 2004), 
surface brightness in i (Ferguson et al. 2007) and GALEX FUV (Munoz-Mateos 
et al. 2007) passbands, and metallicity (Urbaneja et al. 2005; U et al. 2009). 

The gas surface density radial profile derived from observations constrains 
strongly the model parameter space due to large variation of the model gas density 
radial profiles, especially in the inner parts of the galaxy (Fig. 2, panel a). 

The surface brightness radial profile in the i passband constrains well the stellar 
mass of the models. However, we set the same total mass for all models, thus there 
are no signihcant variations in the i passband radial profiles, except in the outer 
parts of the galaxy, where model radial profiles span a wide range (Fig. 2, panel b). 

The GALEX FUV passband surface-brightness radial profile constrains well 
the SF rate along the model galaxy radius (Fig. 2, panel c). However, this profile is 
not very sensitive to the model parameters since the accretion time-scale exceeds 
the gas consumption time-scale, i.e., SF is gas-accretion regulated (Elmegreen 2014 
and references therein). 
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Fig. 3. Same as Fig. 2, but only models with r.m.s. deviations from the observed 
radial profiles (gas surface density and i passband) smaller than 10% and in the radial 
distance range 2-7 kpc are shown. 



Fig. 4. 3-D parameter space of the models shown in Fig. 3. The gray-shaded circles 
indicate the values of a parameter. The small black points at each node of the model 
grid indicate a fully explored extent of the parameter space. 
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X. kpc X. kpc 

Fig. 5. The galaxy models from the “degeneracy valley” (Fig. 4) in the GALEX FUV 
passband: (a) Pt = 0.30, e = 1%; (b) Pt = 0.31, e = 0.6%; (c) Pt = 0.32, e = 0.3%; (d) 
Pt = 0.34, e = 0.2%. For all models, a = 2. The models possess very similar 1-D radial 
profiles, but they demonstrate different 2-D patterns of the young SF regions. 


The metallicity radial profiles (Fig. 2, panel d) are sensitive to the model param¬ 
eters, however, they are of limited use because of large scatter in the observational 
data. 

Therefore, to further constrain the model parameter space we used only the 
gas surface density and the i passband radial profiles. For each model we calcu¬ 
lated the r.m.s. deviations from both observed radial profiles in the most reliable 
(for models) radial distance range, i.e., 2-7 kpc. The inner region of the model 
galaxy has an increasingly anisotropic grid and the outer parts are affected by the 
uncertainty in the critical gas surface density, Ec- In Fig. 3 we show the models 
selected by r.m.s. deviations from the observed radial profiles (gas surface density 
and i passband) being less than 10%. 

The parameter space of the selected models is shown in Fig. 4. A prominent 
feature of this figure is the “degeneracy valley” - the region in the parameter space 
where reside the models selected by the 1-D radial profile procedure described 
above. The “degeneracy valley” is seen because the 1-D radial profiles are derived 
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by convolving the 2-D distributions of SF regions, e.g., SF in the model galaxy 
being smooth and inefficient or patchy and efficient would produce similar 1-D 
radial profiles. 

This effect is shown in Fig. 5, where four different models taken from the 
“degeneracy valley” of Fig. 4 are plotted. Fig. 5 illustrates that, by increasing Pt 
and decreasing e values, SF becomes more smooth. It is challenging to break those 
degeneracies by using only 1-D observed radial profiles, however, feature-rich 2-D 
model structures (Fig. 5) imply the possibility to solve this problem by employing 
spatially resolved observations of real galaxy disks. 

5. CONCLUSIONS 

We presented a 2-D model of galaxy disk evolution, based on the stochastic 
self-propagating star formation process. The parameterized gas and star dynamics, 
together with self-consistent chemical and photometric evolution of stellar popu¬ 
lations, are implemented in the model. The galaxy build-up prescription is taken 
from the N-body DM simulations by Fakhouri et al. (2010). 

We applied the model for the interpretation of observational data on the M33 
galaxy and successfully reproduced the 1-D radial profiles of the gas surface den¬ 
sity, GALEX FUV and i passband surface brightness, and stellar metallicity. We 
explored possible degeneracies of the model parameters by computing the exten¬ 
sive grid of models. We have found that observational data, averaged as 1-D 
radial profiles, cannot fully characterize galaxy evolution. In order to break the 
degeneracies, the analysis of 2-D SF region patterns is additionally required. 
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